#ifndef IMAGE_HPP_
#define IMAGE_HPP_

namespace quadil{

/**
 * Generic struct describing a single chanel digital image.
 * Image is the most essential part in the library.
 * Key consepts:
 * 1)No boundary checking
 * 2)Intended for standard numerical types
 * 3)Intended only for single chanel images (grayscale). Multichanel images can be obtained by hacks
 * 4)Indexing starts from the top left corner of the image
 * 5)Images and Regions Of Interest are implemented as the same class with diferent constructors
 * 6)Very minimal error checking in general
 * 7)An image is imutable
 * 8)An image should have minimum height of 1 and a minimum width of 1
 * 9)For all methods to work as intented, type T must be able to be assigned the value 0.
 * 10)For all methods to work as intented, type T must be able to be to and from stantdard types.
 * 11)All pixels set to 0 are considered background pixels, all other are considered foreground.
 */
template <typename T>struct Image{

	///the number of by
	static const int __memory_allignment__=8;
	///Smart pointer handling memory allocation shold not be handled
	boost::shared_array<T> data;

	///Pointer to the first pixel (top left)
	T* const fp;

	///The width of the image in pixels
	const int width;

	///The height of the image in Pixels
	const int height;

	///The number of sizeof(T) we have to add to go from Pixel(x,y) to Pixel(x,y+1). In normal Images it should be the same as width
	const int lineStride;

	static T* __getFp__(const boost::shared_array<T>& data){
		T* dataptr=data.get();
		return ((T*)((long(dataptr))+(__memory_allignment__-(long(dataptr)%__memory_allignment__))));
	}
	static int __getLineStride__(int width){
		return ((((width*sizeof(T))/__memory_allignment__)+1)*__memory_allignment__)/sizeof(T);
	}
	/**
	 * Default constructor
	 * Allocates memory for the image.
	 * this should be the only constructor allocating
	 * memory.
	 * \param w the desired image width
	 * \param h the desired image height
	 */
	Image(int w,int h):width(w),height(h),
		lineStride(__getLineStride__(w)),
		data(new T[__getLineStride__(w)*h+__memory_allignment__]),
		fp(__getFp__(data))	{}

	/**
	 * Shalow copy constructor
	 * very light weight
	 * new and old image share same data
	 * \param i the original image
	 * */
	Image(const Image<T>& i):width(i.width),height(i.height),lineStride(i.lineStride),data(i.data),fp(i.fp){}

	/**
	 * Unsafe Constructor used to create subimages.
	 * \param i the ParentImage
	 * \param w the desired width of the subimage
	 * \param h the desired height of the subimage
	 * \param firstPixel pointer to the top left pixel
	 * Its up to the caller to make shure w,h and firstPixel are valid:
	 * 1)(i.fp-firstPixel)%i.lineStride+w <= i.width
	 * 2)(i.fp-firstPixel)/i.lineStride+w <= i.height
	 **/
	Image(const Image<T>& i,int w,int h,T *const  firstPixel):data(i.data),width(w),height(h),fp(firstPixel),lineStride(i.lineStride){}

	/**
	 * Creates an image that has a safe border so that accessing pixels outside the image boundaries wont produce errors.
	 * Once the image is created there is no way of knowing the size of the safe region around the image so its up to the
	 * caller to keep track of that region. The rational for using such an image is that some times it is faster and simpler
	 * to iterate outside (but still near) the images edges than having detailed boundary checking. This kind of image is used
	 * \param w the width of the resulting image
	 * \param h the height of the resulting image
	 * \param leftBorder the number of columns beyond the left border (column 0) one can access safely
	 * \param topBorder the number of rows beyond the top border (row 0) one can access safely
	 * \param rightBorder the number of columns beyond the right border (column (width-1)) one can access safely
	 * \param bottomBorder the number of rows beyond the bootom border (row (height-1)) one can access safely
	 * \param val the value to witch all pixels in the container will be set. In most cases 0 is the most convinient
	 */
	static Image<T> createSafeBorderImage(int w,int h,int leftBorder,int topBorder,int rightBorder,int bottomBorder,T val){
		Image<T> container(w+leftBorder+rightBorder,h+topBorder+bottomBorder);
		container.setAll(val);
		/*std::cerr<<"createSafeBorderImage=\n\twidth="<<container.width<<"\n\theight="<<container.height<<"\n\tleft="
				<<leftBorder<<"\n\ttop="<<topBorder<<"\n\tright="<<rightBorder<<"\n\tbottom="<<bottomBorder<<"\n\tgetSubImage("<<leftBorder<<","<<topBorder<<
				","<<leftBorder+container.width-1<<","<<topBorder+container.height-1<<")\n";*/
		return container.getSubImage(leftBorder,topBorder,leftBorder+w-1,topBorder+h-1);
	}
	/**
	 * Creates an image that has a safe border so that accessing pixels outside the image boundaries wont produce errors.
	 * Once the image is created there is no way of knowing the size of the safe region around the image so its up to the
	 * caller to keep track of that region. The rational for using such an image is that some times it is faster and simpler
	 * to iterate outside (but still near) the images edges than having detailed boundary checking. This kind of image is used
	 * \param i the image wicth will be copied in to the "visible" area of the new image
	 * \param leftBorder the number of columns beyond the left border (column 0) one can access safely
	 * \param topBorder the number of rows beyond the top border (row 0) one can access safely
	 * \param rightBorder the number of columns beyond the right border (column (width-1)) one can access safely
	 * \param bottomBorder the number of rows beyond the bootom border (row (height-1)) one can access safely
	 * \param val the value to witch all safe border pixels will be set. In most cases 0 is the most convinient
	 */
	static Image<T> createSafeBorderCopy(Image<T> &i,int leftBorder,int topBorder,int rightBorder,int bottomBorder,T val){
		Image<T> container(i.width+leftBorder+rightBorder,i.height+topBorder+bottomBorder);
		int y,x,lineStrideDiff;
		lineStrideDiff=leftBorder+rightBorder;
		T* curRowOut;
		T* curRowOutStop;
		const T* curRowIn;
		curRowOut=container.getRow(0);
		curRowOutStop=curRowOut*container.lineStride*topBorder;
		while(curRowOut!=curRowOutStop){//initialise top border
			*curRowOut=val;
			curRowOut++;
		}
		for(y=0;y<i.height;y++){
			curRowOutStop=container.getRow(y+topBorder)+leftBorder;
			curRowIn=i.getRow(y);
			while(curRowOut!=curRowOutStop){//initialise left border
				*curRowOut=val;
				curRowOut++;
			}
			for(x=0;x<i.width;x++){//copy image
				curRowOut[x]=curRowIn[x];
			}
			curRowOutStop=container.getRow(y+topBorder+1);
			while(curRowOut!=curRowOutStop){//initialise right border
				*curRowOut=val;
				curRowOut++;
			}
		}
		curRowOutStop=container.getRow(container.height);
		while(curRowOut!=curRowOutStop){//initialise bootom border
			*curRowOut=val;
			curRowOut++;
		}
		return container.getSubImage(leftBorder,topBorder,leftBorder+container.width-1,topBorder+container.height-1);
	}


	/**
	 * Creates an image that contains a rectangular region of the image
	 * Memory is shared beetween *this and the result
	 * \param l the left most position of pixels of the new image
	 * \param t the top most position of pixels of the new image.
	 * \param r the right most position of pixels of the new image.
	 * \param b the bottom most position of pixels of the new image.
	 * \return An image that contains a part the image
	 * for the routine to execute successfully there must:
	 * 0<=l<r<this->width and 0<=t<b<this->height
	 */
	Image<T> getSubImage(unsigned int l,unsigned int t,unsigned int r,unsigned int b){
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
		if(r<=l || b<=t || width<=r || height<=b){
			std::cerr<<"Image::getSubImage("<<l<<","<<t<<","<<r<<","<<b<<") ,height="<<height<<" , width="<<width<<"\n";
			throw "Exception: quadil::image::getSubImage : size missmatch. the sub image rectangle exeeded allowed dimentions";
		}
#endif
		return Image<T>(*this,r-l+1,b-t+1,fp+l+t*lineStride);
	}
	/**
	 * Unsafe Read-Write access to an images row as a pointer to The first pixel of the row
	 * \param rowNum the desired row valid range [0,height)
	 * \return a pointer to the first pixel of the desired row
	**/
	T* const getRow(int rowNum){
		return fp+lineStride*rowNum;
	}
	/**
	 * Unsafe Read Only access to an images row as a pointer to The first pixel of the row
	 * \param rowNum the desired row valid range [0,height)
	 * \return a const pointer to the first pixel of the desired row
	**/
	T const * const getRow(int rowNum)const{
		return fp+lineStride*rowNum;
	}

	/**
	 * Unsafe Read Write access to an images pixel as a reference
	 * For iterations over all pixels of the image image it adds a
	 * significant overhead compaired to the method demonstrated in the
	 * implementation of  Image::binarize
	 * \param x the desired column, valid range [0,width)
	 * \param y the desired row, valid range [0,height)
	 * \return a reference to the desired pixel
	**/
	T& getPixel(int x,int y){
		return *(fp+lineStride*y+x);
	}
	/**
	 * Unsafe Read Only access to an images pixel as a reference
	 * For iterations over all pixels of the image image it adds a
	 * significant overhead compaired to the method menstrated in the
	 * implementation of  Image::binarize
	 * \param x the desired column, valid range [0,width)
	 * \param y the desired row, valid range [0,height)
	 * \return a const reference to the desired pixel
	**/
	const T& getPixel(int x,int y)const{
		return *(fp+lineStride*y+x);
	}

	/**
	 * Sets all pixels to a given value.
	 * The source code of this function is a good example of a simple
	 * iteration over all pixels.
	 * \param val the value all Pixels will be set to
	**/
	void setAll(T val){
		T* curPixel;
		T* curRowEnd;
		int rowNum;
		for(rowNum=0;rowNum<height;rowNum++){
			curPixel=this->getRow(rowNum);
			curRowEnd=curPixel+width;
			while(curPixel!=curRowEnd){
				*curPixel=val;
				curPixel++;
			}
		}
	}
	/**
	 * Binarizes the image.
	 * All pixels bellow the threshold are considered false,
	 * All pixels equal to or grater than the threshold are considered true.
	 * \param minTrueValue the binarization threshold
	 * \param trueValue : all true pixels will be set to trueVal
	 * \param falseValue: all false pixels will be set to falseVal
	 **/
	void binarize(T minTrueValue,T trueVal,T falseVal){
		T* curPixel;
		T* curRowEnd;
		int rowNum;
		for(rowNum=0;rowNum<height;rowNum++){
			curPixel=this->getRow(rowNum);
			curRowEnd=curPixel+width;
			while(curPixel!=curRowEnd){
				*curPixel=(*curPixel<minTrueValue)?falseVal:trueVal;
				curPixel++;
			}
		}
	}
	/**
	 * Creates a new Image with its own storage that contains a copy
	 * of the image. The method is template so the method can be used
	 * as rough conversion between diferent image types.
	 **/
	template <typename K> Image<K> deepCopyGeneric()const{
		Image<K> res(width,height);
		T* const curRowIn;
		K* curRowOut;
		int rowNum,columnNum;
		for(rowNum=0;rowNum<height;rowNum++){
			curRowIn=this->getRow(rowNum);
			curRowOut=res.getRow(rowNum);
			columnNum=0;
			while(columnNum<width){
				curRowOut[columnNum]=curRowIn[columnNum];
				columnNum++;
			}
		}
		return res;
	}
	/**
	 * Creates a new Image with its own storage that contains a copy
	 * of the image.
	 **/
	Image<T> deepCopy()const{
		Image<T> res(width,height);
		const T*  curRowIn;
		T* curRowOut;
		int rowNum,columnNum;
		for(rowNum=0;rowNum<height;rowNum++){
			curRowIn=this->getRow(rowNum);
			curRowOut=res.getRow(rowNum);
			columnNum=0;
			while(columnNum<width){
				curRowOut[columnNum]=curRowIn[columnNum];
				columnNum++;
			}
		}
		return res;
	}
	/**
	 * All pixels in the image that are background in the mask image will be set to background.
	 * the method is tmplated so that any image can be used as mask
	 * \param mask the image witch is used a mask
	 * *this and mask should have the same size.
	 */
	template <typename K> void applyMaskGeneric(const Image<K>& mask){
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
		if(height!=mask.height ||width!=mask.width){
			throw "Exception: quadil::image::applyMaskGeneric : size missmatch. *this and mask should have the same size";
		}
#endif
		const K*  curRowIn;
		T* curRow;
		int rowNum,columnNum;
		for(rowNum=0;rowNum<height;rowNum++){
			curRow=this->getRow(rowNum);
			curRowIn=mask.getRow(rowNum);
			columnNum=0;
			while(columnNum<width){
				(curRowIn[columnNum])||(curRow[columnNum]=0);
				//curRow[columnNum]=curRowIn[columnNum];
				columnNum++;
			}
		}
	}
	/**
	 * Generalization of image inversion
	 * In most cases an image is encoded in a ranges of values
	 * [0 1] or [0 255].
	 * In any case where an image is encoded as a range [0 maxval]
	 * the image is inverted by applying to each pixel the transform newValue=maxval-oldValue
	 * this allows the user of the method to treat the image as an image of any possible encoding
	 * \param the value from witch all pixels will be subtracted. should be equal or greater than maximum value in the image
	 */
	void subtractFrom(T val){
		T* curPixel;
		T* curRowEnd;
		int rowNum;
		for(rowNum=0;rowNum<height;rowNum++){
			curPixel=this->getRow(rowNum);
			curRowEnd=curPixel+width;
			while(curPixel!=curRowEnd){
				*curPixel=val-*curPixel;
				curPixel++;
			}
		}
	}
	/**
	 *Iterates the Image to find the minimum and maximum in the image.
	 * \return an std::pair<T,T> containing minimum and maximum respectively.
	 */
	void getMinMax(T& min,T& max){
		min=*fp;
		max=*fp;
		T* curPixel;
		T* curRowEnd;
		int rowNum;
		for(rowNum=0;rowNum<height;rowNum++){
			curPixel=this->getRow(rowNum);
			curRowEnd=curPixel+width;
			while(curPixel!=curRowEnd){
				(*curPixel>max)&&(max=*curPixel);
				(*curPixel<min)&&(min=*curPixel);
				curPixel++;
			}
		}
	}
	/**
	 * Normalises the image linearly.
	 * It maps values of an image from a value range to an new value range
	 * \param newMin the minimum value in the new value range
	 * \param newMax the maximum value in the new value range
	 * \param oldMin the presumably minimum value in the old value range
	 * \param oldMax the presumably maximum value in the new value range
	 */
	void normaliseLinearly(T newMin,T newMax,T oldMin, T oldMax ){
		double newRange=newMax-newMin;
		double oldRange=oldMax-oldMin;
		T* curPixel;
		T* curRowEnd;
		int rowNum;
		for(rowNum=0;rowNum<height;rowNum++){
			curPixel=this->getRow();
			curRowEnd=curPixel+width;
			while(curPixel!=curRowEnd){
				*curPixel=((*curPixel-oldMin)/oldRange)*newRange+newMin;
				curPixel++;
			}
		}
	}
	/**
	 * Normalises the image linearly.
	 * It maps values of an image from the existing range to an new value range
	 * The operation consists of two iterations:
	 * One iteration to find the existing range
	 * and one iteration to map the calculated range to the new range.
	 * \param newMin the minimum value in the new value range
	 * \param newMax the maximum value in the new value range
	 * \param oldMin the presumably minimum value in the old value range
	 * \param oldMax the presumably maximum value in the new value range
	 */
	void normaliseLinearly(T newMin,T newMax){

		T oldMin,oldMax;
		this->getMinMax(oldMin,oldMax);
		double newRange=newMax-newMin;
		double oldRange=oldMax-oldMin;
		T* curPixel;
		T* curRowEnd;
		int rowNum;
		for(rowNum=0;rowNum<height;rowNum++){
			curPixel=this->getRow(rowNum);
			curRowEnd=curPixel+width;
			while(curPixel!=curRowEnd){
				*curPixel=((*curPixel-oldMin)/oldRange)*newRange+oldMin;
				curPixel++;
			}
		}
	}
	/**
	 * Smears the image horizontaly.
	 * The smearing usualy applies to binary images.
	 * This is a generalisation of smearing to any type, consistent with smearing in binary images.
	 * In this context, Pixels can be viewed one of three classes:
	 * 1)Background pixels (pixels equal to 0)
	 * 2)True pixels (pixels equal to trueVal)
	 * 3)False pixels (non zero pixels that have other values than trueVal)
	 * Smearing will fill continius blocks of background pixels that have True pixels on both their
	 * edges if the size of the continius block is equal or less to maxGapToFill.
	 * example: let 'B' be background pixels, 'T' be true pixels, 'F' be false pixels and maxGapToFill is 5.
	 * a line 'BBTTTBBTBFBTTBBBBBTBBBBBBTTTBBBT' will be changed to
	 *   line 'BBTTTTTTBFBTTTTTTTTBBBBBBTTTTTTT'
	 *
	 * \param maxGapToFill the maximum size of continius background pixels that will be filled
	 * \param trueVal the value that filled background pixels will be set to.
	 */
	void applyHorizontalRLSA(int maxGapToFill,T trueVal){
		T* curLine;
		int x,y,i;
		int lastTruePixel;
		for(y=0;y<height;y++){
			curLine=fp+y*lineStride;
			lastTruePixel=(curLine[0]==trueVal)?0:-1;
			x=1;
			while(x<width){
				if(curLine[x]!=0){
					if(curLine[x]!=trueVal){
						lastTruePixel=-1;
					}else{

						if(curLine[x-1]!=trueVal && lastTruePixel>=0 && (x-lastTruePixel)<=maxGapToFill){
							for(i=x-1;i>lastTruePixel;i--){

								curLine[i]=trueVal;
							}
						}
						lastTruePixel=x;
					}
				}
				x++;
			}
		}
	}
	/**
	 * Smears the image verticaly.
	 * for more info look at Image::applyHorizontalRLSA
	 * \param maxGapToFill the maximum size of continius background pixels that will be filled
	 * \param trueVal the value that filled background pixels will be set to.
	 */
	void applyVerticalRLSA(int maxGapToFill,T trueVal){
		T* curColumn;
		int x,y,i,added=0;
		int lastTruePixel;
		for(x=0;x<width;x++){
			curColumn=fp+x;
			lastTruePixel=(curColumn[0]==trueVal)?0:-1;
			y=0;
			while(y<height*lineStride){
				if(curColumn[y]!=0){
					if(curColumn[y]!=trueVal){
						lastTruePixel=-1;
					}else{
						if(curColumn[y-lineStride]!=trueVal && lastTruePixel>=0 && (y-lastTruePixel)<=maxGapToFill*lineStride){
							for(i=y-lineStride;i>lastTruePixel;i-=lineStride){
								curColumn[i]=trueVal;
								added++;
							}
						}
						lastTruePixel=y;
					}
				}
				y+=lineStride;
			}
		}
	}
	/**
	 * Smears the image diagonaly along the (Right Top)-(Left Bottom) diagonal.
	 * for more info look at Image::applyHorizontalRLSA
	 * \param maxGapToFill the maximum size of continius background pixels that will be filled
	 * \param trueVal the value that filled background pixels will be set to.
	 */
	void applyDiagonalRTLB_RLSA(int maxGapToFill,T trueVal){
		T* curColumn;
		int x,y,tmpx,tmpy,i;
		int lastTruePixel;
		for(x=0;x<width;x++){
			curColumn=fp+x;
			lastTruePixel=-1;
			y=0;
			tmpx=x;
			while((y<height*lineStride)&&(tmpx>0) ){
				if(curColumn[y]!=0){
					//std:cerr<<"([y="<<y/lineStride;
					if(curColumn[y]!=trueVal){
						lastTruePixel=-1;
					}else{
						if(curColumn[y-(lineStride-1)]!=trueVal && lastTruePixel>=0 && (y-lastTruePixel)<=maxGapToFill*lineStride){
							for(i=y;i>lastTruePixel;i-=(lineStride-1)){
								curColumn[i]=trueVal;
							}
						}
						lastTruePixel=y;
					}
				}
				y+=(lineStride-1);
				tmpx--;
			}
		}
		for(tmpy=0;tmpy<height;tmpy++){
			curColumn=fp+width-1;
			lastTruePixel=-1;
			y=tmpy*lineStride+width-1;
			tmpx=width-1;
			while((y<height*lineStride)&&(tmpx>0) ){
				if(curColumn[y]!=0){
					if(curColumn[y]!=trueVal){
						lastTruePixel=-1;
					}else{
						if(curColumn[y-(lineStride-1)]!=trueVal && lastTruePixel>=0 && (y-lastTruePixel)<=maxGapToFill*lineStride){
							for(i=y;i>lastTruePixel;i-=(lineStride-1)){
								curColumn[i]=trueVal;
							}
						}
						lastTruePixel=y;
					}
				}
				y+=(lineStride-1);
				tmpx--;
			}
		}
	}
	/**
	 * Smears the image diagonaly along the (Left Top)-(Right Bottom) diagonal.
	 * for more info look at Image::applyHorizontalRLSA
	 * \param maxGapToFill the maximum size of continius background pixels that will be filled
	 * \param trueVal the value that filled background pixels will be set to.
	 */
	void applyDiagonalLTRB_RLSA(int maxGapToFill,T trueVal){
		T* curColumn;
		int x,y,tmpx,tmpy,i;
		int lastTruePixel;
		for(x=0;x<width;x++){
			curColumn=fp+x;
			lastTruePixel=-1;
			y=0;
			tmpx=x;
			while((y<height*lineStride)&&(tmpx<width-1) ){
				if(curColumn[y]!=0){
					if(curColumn[y]!=trueVal){
						lastTruePixel=-1;
					}else{
						if(curColumn[y-(lineStride+1)]!=trueVal && lastTruePixel>=0 && (y-lastTruePixel)<=maxGapToFill*lineStride){
							for(i=y;i>lastTruePixel;i-=(lineStride+1)){
								curColumn[i]=trueVal;
							}
						}
						lastTruePixel=y;
					}
				}
				y+=(lineStride+1);
				tmpx++;
			}
		}
		for(tmpy=1;tmpy<height;tmpy++){
			curColumn=fp;
			y=tmpy*lineStride;
			tmpx=0;
			lastTruePixel=-1;
			while((y<height*lineStride)&&(tmpx<width-1) ){
				if(curColumn[y]!=0){
					if(curColumn[y]!=trueVal){
						lastTruePixel=-1;
					}else{
						if(curColumn[y-(lineStride+1)]!=trueVal && lastTruePixel>=0 && (y-lastTruePixel)<=maxGapToFill*lineStride){
							for(i=y;i>lastTruePixel;i-=(lineStride+1)){
								curColumn[i]=trueVal;
							}
						}
						lastTruePixel=y;
					}
				}
				y+=(lineStride+1);
				tmpx++;
			}
		}
	}
	/**
	 * Produces an Image that contains all connected components labeled as diferent values
	 * This is a from scratch implementation that is implemented with 2 iterations of the image.
	 * A connected component is considered to be any set of foreground pixels that touch in the
	 * 8 neighbors neighborhood.
	 * \return an Image<t_label> having the same size as *this.
	 */
	void labelConnectedComponents8(Image<t_label>&out)const{
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
		if(out.width!=width ||out.height!=height){
			throw "Exception: quadil::image::labelConnectedComponents8 : size missmatch. *this and out should have the same size";
		}
#endif
		Image<T> i(*this);
		__equivalence__::EquivalenceSets equiv;
		int inSum;
		T* inNorth;
		T* inWest;
		T* inLine;
		T* inNorthWest;
		T* inNorthEast;
		t_label* outNorth;
		t_label* outWest;
		t_label* outLine;
		t_label* outNorthWest;
		t_label* outNorthEast;
		int x,y,inLineStride,outLineStride;
		inLineStride=i.lineStride;
		outLineStride=out.lineStride;
		inLine=(T*)(i.getRow(0));
		inWest=inLine-1;
		outLine=out.getRow(0);
		outWest=outLine-1;
		if(inLine[0]){
			outLine[0]=equiv.getNewCode();
		}else{
			outLine[0]=0;
		}
		for(x=1;x<width;x++){
			if(inLine[x]){
				if(inWest[x]){
					outLine[x]=outWest[x];
				}else{
					outLine[x]=equiv.getNewCode();
				}
			}else{
				outLine[x]=0;
			}
		}
		for(y=1;y<height;y++){
			inLine=(T*)i.getRow(y);
			inWest=inLine-1;
			inNorth=inLine-i.lineStride;
			inNorthEast=inNorth+1;
			inNorthWest=inNorth-1;
			outLine=out.getRow(y);
			outWest=outLine-1;
			outNorth=outLine-(out.lineStride);
			outNorthEast=outNorth+1;
			outNorthWest=outNorth-1;
			if(inLine[0]){
				if(inNorth[0]||inNorthEast[0]){
					if(inNorth[0]&& inNorthEast[0]){
						if(outNorth[0]!=outNorthEast[0]){
							equiv.setEquivalent(outNorthEast[0],outNorth[0]);
						}
						outLine[0]=outNorth[0];
					}else{
						outLine[0]=outNorth[0]+outNorthEast[0];
					}
				}else{
					outLine[0]=equiv.getNewCode();
				}
			}
			for(x=1;x<width-1;x++){
				if(inLine[x]){
					inSum=(inWest[x]!=0)+(inNorthWest[x]!=0)*2+(inNorth[x]!=0)*4+(inNorthEast[x]!=0)*8;
					switch (inSum){
						case 15:
							if(outWest[x]!=outNorth[x]){
								equiv.setEquivalent(outWest[x],outNorth[x]);
							}
							if(outWest[x]!=outNorthWest[x]){
								equiv.setEquivalent(outWest[x],outNorthWest[x]);
							}
							if(outWest[x]!=outNorthEast[x]){
								equiv.setEquivalent(outWest[x],outNorthEast[x]);
							}
							outLine[x]=outWest[x];break;
						case 0:
							outLine[x]=equiv.getNewCode();break;
						case 1:
						case 2:
						case 4:
						case 8:
							outLine[x]=outWest[x]+outNorth[x]+outNorthWest[x]+outNorthEast[x];break;
						case 3:
							if(outWest[x]!=outNorthWest[x]){
								equiv.setEquivalent(outWest[x],outNorthWest[x]);
							}
							outLine[x]=outWest[x];break;
						case 6:
							if(outNorth[x]!=outNorthWest[x]){
								equiv.setEquivalent(outNorth[x],outNorthWest[x]);
							}
							outLine[x]=outNorthWest[x];break;
						case 12:
							if(outNorth[x]!=outNorthEast[x]){
								equiv.setEquivalent(outNorth[x],outNorthEast[x]);
							}
							outLine[x]=outNorth[x];break;
						case 5:
							if(outNorth[x]!=outWest[x]){
								equiv.setEquivalent(outWest[x],outNorth[x]);
							}
							outLine[x]=outWest[x];break;
						case 10:
							if(outNorthWest[x]!=outNorthEast[x]){
								equiv.setEquivalent(outNorthEast[x],outNorthWest[x]);
							}
							outLine[x]=outNorthWest[x];break;
						case 9:
							if(outWest[x]!=outNorthEast[x]){
								equiv.setEquivalent(outNorthEast[x],outWest[x]);
							}
							outLine[x]=outWest[x];break;
						case 7:
							if(outWest[x]!=outNorthWest[x]){
								equiv.setEquivalent(outNorthWest[x],outWest[x]);
							}
							if(outWest[x]!=outNorth[x]){
								equiv.setEquivalent(outNorth[x],outWest[x]);
							}
							outLine[x]=outWest[x];break;
						case 13:
							if(outWest[x]!=outNorth[x]){
								equiv.setEquivalent(outNorth[x],outWest[x]);
							}
							if(outWest[x]!=outNorthEast[x]){
								equiv.setEquivalent(outNorthEast[x],outWest[x]);
							}
							outLine[x]=outWest[x];break;
						case 11:
							if(outWest[x]!=outNorthWest[x]){
								equiv.setEquivalent(outNorthWest[x],outWest[x]);
							}
							if(outWest[x]!=outNorthEast[x]){
								equiv.setEquivalent(outNorthEast[x],outWest[x]);
							}
							outLine[x]=outWest[x];break;
						case 14:
							if(outNorth[x]!=outNorthWest[x]){
								equiv.setEquivalent(outNorthWest[x],outNorth[x]);
							}
							if(outNorth[x]!=outNorthEast[x]){
								equiv.setEquivalent(outNorthEast[x],outNorth[x]);
							}
							outLine[x]=outNorthWest[x];break;
					}
				}else{
					outLine[x]=0;
				}
			}
			if(inLine[x]){
				inSum=(inWest[x]!=0)+(inNorthWest[x]!=0)*2+(inNorth[x]!=0)*4;
				switch (inSum){
					case 0:
						outLine[x]=equiv.getNewCode();break;
					case 1:
					case 2:
					case 4:
						outLine[x]=outWest[x]+outNorth[x]+outNorthWest[x];break;
					case 3:
						if(outWest[x]!=outNorthWest[x]){
							equiv.setEquivalent(outWest[x],outNorthWest[x]);
						}
						outLine[x]=outWest[x];break;
					case 6:
						if(outNorth[x]!=outNorthWest[x]){
							equiv.setEquivalent(outNorth[x],outNorthWest[x]);
						}
						outLine[x]=outNorthWest[x];break;
					case 5:
						if(outNorth[x]!=outWest[x]){
							equiv.setEquivalent(outWest[x],outNorth[x]);
						}
						outLine[x]=outWest[x];break;
					case 7:
						if(outWest[x]!=outNorthWest[x]){
							equiv.setEquivalent(outNorthWest[x],outWest[x]);
						}
						if(outWest[x]!=outNorth[x]){
							equiv.setEquivalent(outNorth[x],outWest[x]);
						}
						outLine[x]=outWest[x];break;
					}
				}else{
					outLine[x]=0;
				}
		}
		equiv.finalize();
		for(y=1;y<height;y++){
			outLine=out.getRow(y);
			for(x=0;x<width;x++){
				(outLine[x])&&(outLine[x]=equiv[outLine[x]]);
			}
		}
	}
	/**
	 * Produces an Image that contains all connected components labeled as diferent values
	 * This is a from scratch implementation that is implemented with 2 iterations of the image.
	 * A connected component is considered to be any set of foreground pixels that touch in the
	 * 8 neighbors neighborhood.
	 * \return an Image<t_label> having the same size as *this.
	 */
	void labelConnectedComponents4(Image<t_label>& out)const{
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
		if(out.width!=width ||out.height!=height){
			throw "Exception: quadil::image::labelConnectedComponents4 : size missmatch. *this and out should have the same size";
		}
#endif
		Image<T> i(*this);
		__equivalence__::EquivalenceSets equiv;
		T*	inNorth;
		T* inWest;
		T* inLine;
		t_label* outNorth;
		t_label* outWest;
		t_label* outLine;
		int x,y,inLineStride,outLineStride;
		inLineStride=i.lineStride;
		outLineStride=out.lineStride;
		inLine=i.getRow(0);
		inWest=inLine-1;
		outLine=(out.getRow(0));
		outWest=outLine-1;
		if(inLine[0]){
			outLine[0]=equiv.getNewCode();
		}else{
			outLine[0]=0;
		}
		for(x=1;x<width;x++){
			if(inLine[x]){
				if(inWest[x]){
					outLine[x]=outWest[x];
				}else{
					outLine[x]=equiv.getNewCode();
				}
			}else{
				outLine[x]=0;
			}
		}
		for(y=1;y<height;y++){
			inLine=i.getRow(y);
			inWest=inLine-1;
			inNorth=inLine-i.lineStride;
			outLine=out.getRow(y);
			outWest=outLine-1;
			outNorth=outLine-(out.lineStride);
			if(inLine[0]){
				if(inNorth[0]){
					outLine[0]=outNorth[0];
				}else{
					outLine[0]=equiv.getNewCode();
				}
			}else{
				outLine[0]=0;
			}
			for(x=1;x<width;x++){
				if(inLine[x]){
					if(inNorth[x] &&(outNorth[x]==outWest[x])){//both neibors true and same;
						outLine[x]=outNorth[x];
					}else{//if neibors differ or both neibors false
						if(inNorth[x]&&inWest[x]){//both neibors true but diferent;
							outLine[x]=outWest[x];
							equiv.setEquivalent(outWest[x],outNorth[x]);
						}else{//one or on neibors true
							if(inNorth[x]||inWest[x]){//one true neibor
								outLine[x]=outWest[x]+outNorth[x];
							}else{//no true neibors
								outLine[x]=equiv.getNewCode();
							}
						}

					}
				}else{
					outLine[x]=0;
				}
			}
		}
		equiv.finalize();
		for(y=1;y<height;y++){
			outLine=(t_label*)out.getRow(y);
			for(x=0;x<width;x++){
				(outLine[x])&&(outLine[x]=equiv[outLine[x]]);
			}
		}
	}

	/**
	 * method that writes the accumulative sum of *this to out.
	 * out should have the same size as this. The method is templated but
	 * the caller should make shure the type of the accumulator is such that
	 * it can represent the sum: eg Image<t_byte> should be accumulated at least by
	 * an Image<int>. An  Image<int> should be accumulated by an Image<long long int> etc.
	 * The accumulation starts from the top left pixel.
	 * Every pixel out[x,y] will be equal sum(this[0-x,0-y]).
	 * \param out the image in witch the the accumulative sum will be stored
	 */
	template <typename K>void getAccumulativeSum(Image<K>& out,int goBeyondX=0,int goBeyondY=0)const{
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
		if(height!=out.height ||width!=out.width){
			throw "Exception: quadil::image::getAccumulativeSum : size missmatch. *this and out should have the same size";
		}
#endif
		K* curRowOut;
		K* curNorthOut;
		const T* curRowIn;
		int x,y,bx,by;
		K rowSum;
		rowSum=0;
		curRowIn=this->getRow(0);
		curRowOut=out.getRow(0);
		for(x=0;x<width;x++){
			rowSum+=curRowIn[x];
			curRowOut[x]=rowSum;
		}
		bx=0;
		while(bx<goBeyondX){
			curRowOut[x+bx]=curRowOut[x-1];
			bx++;
		}
		for(y=1;y<height;y++){
			rowSum=0;
			curRowIn=this->getRow(y);
			curRowOut=out.getRow(y);
			curNorthOut=out.getRow(y-1);
			for(x=0;x<width;x++){
				rowSum+=curRowIn[x];
				curRowOut[x]=curNorthOut[x]+rowSum;
			}
			bx=0;
			while(bx<goBeyondX){
				curRowOut[x+bx]=curRowOut[x-1];
				bx++;
			}
		}
		curNorthOut=curRowOut;
		while(y<height+goBeyondY){
			curRowOut=out.getRow(y);
			x=width+goBeyondX-1;
			while(x>=0){
				curRowOut[x]=curNorthOut[x];
				x--;
			}
			y++;
		}
		for(y=-goBeyondY;y<height+goBeyondY;y++){
			curRowOut=out.getRow(y)-goBeyondX;
			for(x=0;x<goBeyondX;x++){
				curRowOut[x]=0;
			}
		}
		K* curPixelOut=out.getRow(-goBeyondY)-goBeyondX;
		K* pixelStop=out.getRow(0);
		while(curPixelOut!=pixelStop){
			*curPixelOut=0;
			curPixelOut++;
		}
	}
	static const int ABS_DERIVATION_ASCEND= 1;
	static const int ABS_DERIVATION_DESCEND= 2;
	static const int ABS_DERIVATION_BOTH=3;
	void getAbsDerivativeHorizontal(Image<T>& out,int mode=ABS_DERIVATION_BOTH)const{
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
		if(height!=out.height ||width!=out.width){
			throw "Exception: quadil::image::getAbsDerivativeHorizontal : size missmatch. *this and out should have the same size";
		}
#endif
		int x,y;
		const T* inLine;
		const T* inLineShifted;
		T* outLine;
		if(mode==ABS_DERIVATION_BOTH){
			for (y=0;y<height;y++){
				inLine=this->getRow(y);
				inLineShifted=inLine-1;
				outLine=out.getRow(y);
				outLine[0]=0;
				for(x=1;x<width;x++){
					outLine[x]=(inLine[x]<inLineShifted[x])?(inLineShifted[x]-inLine[x]):(inLine[x]-inLineShifted[x]);
				}
			}
		}
		if(mode==ABS_DERIVATION_ASCEND){
			for (y=0;y<height;y++){
				inLine=this->getRow(y);
				inLineShifted=inLine-1;
				outLine=out.getRow(y);
				outLine[0]=0;
				for(x=1;x<width;x++){
					outLine[x]=(inLine[x]<inLineShifted[x])?(0):(inLine[x]-inLineShifted[x]);
				}
			}
		}
		if(mode==ABS_DERIVATION_DESCEND){
			for (y=0;y<height;y++){
				inLine=this->getRow(y);
				inLineShifted=inLine-1;
				outLine=out.getRow(y);
				outLine[0]=0;
				for(x=1;x<width;x++){
					outLine[x]=(inLine[x]<inLineShifted[x])?(inLineShifted[x]-inLine[x]):(0);
				}
			}
		}

	}
	void getAbsDerivativeVertical(Image<T>& out,int mode=ABS_DERIVATION_BOTH)const{
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
		if(height!=out.height ||width!=out.width){
			throw "Exception: quadil::image::getAbsDerivativeVertical : size missmatch. *this and out should have the same size";
		}
#endif
		int x,y;
		const T* inLine;
		const T* inPreveusLine;
		T* outLine;
		outLine=out.getRow(0);
		for(x=0;x<width;x++){
			outLine[x]=0;
		}
		if(mode==ABS_DERIVATION_BOTH){
			for (y=1;y<height;y++){
				inLine=this->getRow(y);
				inPreveusLine=this->getRow(y-1);
				outLine=out.getRow(y);
				for(x=0;x<width;x++){
					outLine[x]=(inLine[x]<inPreveusLine[x])?(inPreveusLine[x]-inLine[x]):(inLine[x]-inPreveusLine[x]);
				}
			}
		}
		if(mode==ABS_DERIVATION_ASCEND){
			for (y=1;y<height;y++){
				inLine=this->getRow(y);
				inPreveusLine=this->getRow(y-1);
				outLine=out.getRow(y);
				for(x=0;x<width;x++){
					outLine[x]=(inLine[x]<inPreveusLine[x])?(0):(inLine[x]-inPreveusLine[x]);
				}
			}
		}
		if(mode==ABS_DERIVATION_DESCEND){
			for (y=1;y<height;y++){
				inLine=this->getRow(y);
				inPreveusLine=this->getRow(y-1);
				outLine=out.getRow(y);
				for(x=0;x<width;x++){
					outLine[x]=(inLine[x]<inPreveusLine[x])?(inPreveusLine[x]-inLine[x]):(0);
				}
			}
		}
	}

	void getAbsDerivativeDiagonalTLBR(Image<T>& out,int mode=ABS_DERIVATION_BOTH)const{
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
		if(height!=out.height ||width!=out.width){
			throw "Exception: quadil::image::getAbsDerivativeVertical : size missmatch. *this and out should have the same size";
		}
#endif
		int x,y;
		const T* inLine;
		const T* inPreveusLine;
		T* outLine;
		outLine=out.getRow(0);
		for(x=0;x<width;x++){
			outLine[x]=0;
		}
		if(mode==ABS_DERIVATION_BOTH){
			for (y=1;y<height;y++){
				inLine=this->getRow(y);
				inPreveusLine=this->getRow(y-1)-1;
				outLine=out.getRow(y);
				outLine[0]=0;
				for(x=1;x<width;x++){
					outLine[x]=(inLine[x]<inPreveusLine[x])?(inPreveusLine[x]-inLine[x]):(inLine[x]-inPreveusLine[x]);
				}
			}
		}
		if(mode==ABS_DERIVATION_ASCEND){
			for (y=1;y<height;y++){
				inLine=this->getRow(y);
				inPreveusLine=this->getRow(y-1)-1;
				outLine=out.getRow(y);
				outLine[0]=0;
				for(x=1;x<width;x++){
					outLine[x]=(inLine[x]<inPreveusLine[x])?(0):(inLine[x]-inPreveusLine[x]);
				}
			}
		}
		if(mode==ABS_DERIVATION_DESCEND){
			for (y=1;y<height;y++){
				inLine=this->getRow(y);
				inPreveusLine=this->getRow(y-1)-1;
				outLine=out.getRow(y);
				outLine[0]=0;
				for(x=1;x<width;x++){
					outLine[x]=(inLine[x]<inPreveusLine[x])?(inPreveusLine[x]-inLine[x]):(0);
				}
			}
		}

	}
	void getAbsDerivativeDiagonalTRBL(Image<T>& out,int mode=ABS_DERIVATION_BOTH)const{
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
		if(height!=out.height ||width!=out.width){
			throw "Exception: quadil::image::getAbsDerivativeVertical : size missmatch. *this and out should have the same size";
		}
#endif
		int x,y;
		const T* inLine;
		const T* inPreveusLine;
		T* outLine;
		outLine=out.getRow(0);
		for(x=0;x<width;x++){
			outLine[x]=0;
		}
		if(mode==ABS_DERIVATION_BOTH){
			for (y=1;y<height;y++){
				inLine=this->getRow(y);
				inPreveusLine=this->getRow(y-1)+1;
				outLine=out.getRow(y);
				outLine[0]=0;
				for(x=width-2;x>=0;x--){
					outLine[x]=(inLine[x]<inPreveusLine[x])?(inPreveusLine[x]-inLine[x]):(inLine[x]-inPreveusLine[x]);
				}
			}
		}
		if(mode==ABS_DERIVATION_ASCEND){
			for (y=1;y<height;y++){
				inLine=this->getRow(y);
				inPreveusLine=this->getRow(y-1)+1;
				outLine=out.getRow(y);
				outLine[0]=0;
				for(x=width-2;x>=0;x--){
					outLine[x]=(inLine[x]<inPreveusLine[x])?(0):(inLine[x]-inPreveusLine[x]);
				}
			}
		}
		if(mode==ABS_DERIVATION_DESCEND){
			for (y=1;y<height;y++){
				inLine=this->getRow(y);
				inPreveusLine=this->getRow(y-1)+1;
				outLine=out.getRow(y);
				outLine[0]=0;
				for(x=width-2;x>=0;x--){
					outLine[x]=(inLine[x]<inPreveusLine[x])?(inPreveusLine[x]-inLine[x]):(0);
				}
			}
		}

	}

	template <int NBIMAGES> void accumulateWeightedSum(const Image<T>*in,double w,...){
		const Image<T>* inImgs[NBIMAGES];
		double imgWeights[NBIMAGES];
		const T* inLines[NBIMAGES];
		T* outLine;
		//std::cerr<<"Accum 1\n";
		va_list vl;
		va_start(vl,w);
		//std::cerr<<"Accum 2\n";
		int k,x,y;
		inImgs[0]=in;
		imgWeights[0]=w;
		for (k=1;k<NBIMAGES;k++){
			//std::cerr<<"Accum 2.1\n";
			inImgs[k]=va_arg(vl,const Image<T>*);
			//std::cerr<<"Accum 2.2\n";
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
			if(height!=inImgs[k]->height ||width!=inImgs[k]->width){
				throw "Exception: quadil::image::accumulateWeightedSum : size missmatch. *this and all input images should have the same size";
			}
			//std::cerr<<"Accum 2.5\n";
#endif
			//std::cerr<<"Accum 2.6\n";
			imgWeights[k]=va_arg(vl,double);
			//std::cerr<<"Accum 2.7\n";
		}
		//std::cerr<<"Accum 3\n";
		va_end(vl);
		//std::cerr<<"Accum 4\n";
		for(y=0;y<height;y++){
			for(k=0;k<NBIMAGES;k++){
				inLines[k]=inImgs[k]->getRow(y);
			}
			outLine=this->getRow(y);
			for(x=0;x<width;x++){
				outLine[x]=inLines[0][x]*imgWeights[0];
				for(k=1;k<NBIMAGES;k++){
					outLine[x]+=inLines[k][x]*imgWeights[k];
				}
			}
		}
	}
	template <int NBIMAGES> void accumulateWeightedProd(const Image<T>*in,double w,...){
		const Image<T>* inImgs[NBIMAGES];
		double imgWeights[NBIMAGES];
		const T* inLines[NBIMAGES];
		T* outLine;
		std::cerr<<"Accum 1\n";
		va_list vl;
		va_start(vl,w);
		std::cerr<<"Accum 2\n";
		int k,x,y;
		inImgs[0]=in;
		imgWeights[0]=w;
		for (k=1;k<NBIMAGES;k++){
			std::cerr<<"Accum 2.1\n";
			inImgs[k]=va_arg(vl,const Image<T>*);
			std::cerr<<"Accum 2.2\n";
#ifdef QUADIL_TUNE_CHECK_SAME_SIZE_IMAGES
			if(height!=inImgs[k]->height ||width!=inImgs[k]->width){
				throw "Exception: quadil::image::accumulateWeightedSum : size missmatch. *this and all input images should have the same size";
			}
			std::cerr<<"Accum 2.5\n";
#endif
			std::cerr<<"Accum 2.6\n";
			imgWeights[k]=va_arg(vl,double);
			std::cerr<<"Accum 2.7\n";
		}
		std::cerr<<"Accum 3\n";
		va_end(vl);
		std::cerr<<"Accum 4\n";
		for(y=0;y<height;y++){
			for(k=0;k<NBIMAGES;k++){
				inLines[k]=inImgs[k]->getRow(y);
			}
			outLine=this->getRow(y);
			for(x=0;x<width;x++){
				outLine[x]=inLines[0][x]*imgWeights[0];
				for(k=1;k<NBIMAGES;k++){
					outLine[x]*=inLines[k][x]*imgWeights[k];
				}
			}
		}
	}
};


}
#endif /*IMAGE_HPP_*/
